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We suggest a new approach to the detection of gravitational waves using observa- 
tions of a group of millisecond pulsars. In contrast to the usual method, based on 
increasing the accuracy of the arrival times of pulses by excluding possible distorting 
factors, our method supposes that the additive phase noise that is inevitably present 
even in the most accurate observational data has various spectral components, which 
have characteristic amplitudes and begin to appear on different time scales. We use 
the "Caterpillar" (Singular Spectral Analysis, SSA) method to decompose the signal 
into its components. Our initial data are the residuals of the pulse arrival times for 
six millisecond pulsars. We constructed the angular correlation function for compo- 
nents of the decomposition of a given number, whose theoretical form for the case 
of an isotropic and homogeneous gravitational-wave background is known. The indi- 
vidual decomposition components show a statistically significant agreement with the 
theoretical expectations (correlation coefficient p = 0.92 ± 0.10). 
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1. INTRODUCTION 

Long-term analyses of pulsar-timing data have shown that the fractional instability of the 
rotations of some pulsars are comparable to the instabilities of atomic frequency standards, 
and can reach Au/u ~ 10~ 15 . This makes it possible to apply timing data to various 
astronomical and metrological problems. In the present paper, we describe a method for 
detecting gravitational waves using observations of a group of millisecond pulsars. The 
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method is based on the so-called two-point (angular) correlation function [1-3] 

C(6>) = |stog(s)-| + ± + ^(s), (1) 

where x — (1 — cos#)/2, is the angular separation of the pulsars. After the publication of 
[1] the principle of gravitational-wave detection using pulsar-timing data became clear, and 
observations employing this method were initiated in many observatories. 

The main attention in these programs was focused on increasing the accuracy of pulse 
time-of- arrival (TOA) determinations and excluding all possible distorting physical factors 
(instability of the reference clock, influence of the interstellar medium, the intrinsic activity 
of the pulsars themselves, etc.). However, it is clear that increasing the accuracy of the obser- 
vations is a necessary, but not sufficient, condition, since increasingly finer effects are found, 
which have a stochastic nature and are able to influence the gravitational- wave background. 

We suggest here a fundamentally different approach to detecting gravitational waves. It 
is known [4] that physical objects such as atomic time standards or pulsars display noise 
over a broad frequency range. For example, the phase noise of a frequency standard or a 
pulsar contains both white noise, which is usually identified with errors of the detectors, 
and correlated (red) noise with various spectral indexes [5], which begins to appear in data 
collected on various time intervals. Moreover, the analysis of time series corresponding to 
the different types of noises shows that they have characteristic features that can be used 
to discriminate between them. These features suggest that there should exist a method 
for decomposing time series into components having different spectral compositions. If one 
such component is dominated by the gravitational- wave background, the presence of this 
background should be distinguishable in the angular correlation function. We choose the 
"Caterpillar" or Singular Spectral Analysis (SSA) method [6] to decompose the time series 
into separate components. This method is functionally independent, since the time series 
itself is used to obtain the orthogonal basis functions. 

2. OBSERVATIONS 

Observations of the pulsars PSR J0613-0200, J1640+2224, J1643-1224, J1713+0747, 
J1939+2134, and J2145-0750 were carried out using the fully steerable 64-m radio telescope 
of the Kalyazin Radio Astronomical Observatory of the Astro Space Center of the Lebedev 
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Physical Institute [7-10]. The pulses were accumulated in two circular polarizations using a 
spectrum analyzer with an 80-channel filter bank with a bandwidth of 40kHz per channel in 
each polarization [11]. Each pulsar was observed, on average, once in two weeks. The total 
signal sampling time was close to two hours per session. The pulse TOAs were determined 
using a local time standard with accuracy better than 100 ns. The local time standard was 
tied to UTC(SU) via a television channel and to UTC(USNO) via a GPS receiver. The 
pulses were accumulated in three-minute cycles synchronous with the pulsar rotation, with 
the subsequent restart of the system using the newly computed observed pulsar period. The 
data for each three-minute cycle were recorded in a separate file. The pulse delays due to 
frequency dispersion in the individual channels were compensated for in the subsequent data 
reduction. 

The topocentric pulse TOAs were determined by matching the pulsar profile summed 
over a session to a standard high signal-to- noise profile (template). The barycentric pulse 
TOAs, TOA residuals, and refined astrometric and rotational parameters of the pulsars were 
calculated via a least-squares minimization of the TOA residuals using the Tempo package 
[12]. The DD model [13] was used to refine the orbital parameters of the pulsars. We used 
the values given in the pulsar catalog [14] as the initial values. 

Table 1 lists the least-squares-fit timing parameters for the pulsars: the right ascension 
a(J2000), declination 5(J2000), proper motions fx a , /j, s (in mas/yr) for epoch J2000, the 
rotational frequency / and its derivative / (in s _1 and s -2 ), the dispersion measure DM 
(in pc/cm 3 ), the projection of the semi-major axis of the pulsar orbit onto the line of sight 
x (in light seconds), the eccentricity of the orbit e, the epoch of the pulsar's transit Tn (in 
MJD), the longitude of the orbit periastron u (in deg), the rms error of the residuals after 
improving the pulsar timing parameters a (in /xs). Since the observations were carried out 
in a single-frequency regime, we did not correct for the DM . The formal least-squares error 
of the last significant digit of the parameter is given in parentheses. 

To compute the components of the pulse-TOA decompositions using the SSA method, 
the time-series were averaged over 40-day intervals, filling gaps via a linear interpolation 
between adjacent measurements. We took only common parts of the pulse-TOA series in 
the range MJD = 51000-53000. After averaging and interpolation, the pulse-TOA series 
contained iV = 51 points. The series of pulse-TOA residuals averaged as described above 
are shown in Fig. 1. 
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3. THE "CATERPILLAR" (SSA) METHOD 

We will now describe the steps in the SSA method for a one- dimensional time-series 
{xi}f =1 with length iV following [15]. 



1. The first step is unfolding the one-dimensional series into a multi-dimensional series. 
We take a number M < N (the length of "caterpillar"), specify k = N — M + 1 and 



form the matrix X = (x^) 1 *'^ with elements Xij = x i+ j-\. 



( 
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2. Next, the following matrix is solved: 



R = XX T , 



(3) 



where () T denote conjugation. 

3. The eigenvalues and eigenvectors of the matrix R are then computed, i.e., we find 
expand it in the form 

R = PAP T , 



where 
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is a diagonal matrix of the eigenvalues and 
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is an orthogonal matrix of the eigenvectors of the matrix R. Note that [15]: P T — P 1 , 
pTp = pp T = /mj A = P T Rp ^ ^« \. = M; Y\f =1 \ = detR, where I M is the 

identity matrix with dimension M and det(-) is the determinant of the matrix (■). 

4. The transformation to the principal components is accomplished using the formula 

XP = Y = (y u y 2 ,...,y M ), (7) 
and the inverse transformation to the matrix X is computed X = YP T . 

5. Further is the reconstruction of the initial series. The reconstruction procedure is 
based on the formula X = Y*P T . It is said that the reconstruction is carried out using 
a given set of principal components if the matrix Y* is obtained from the matrix Y by 
setting all components not included in the set of principle components to zero. Thus, 
we can obtain an approximation for the matrix of interest, or to some interpretable part 
of this matrix. The one-dimensional time-series is obtained via a diagonal averaging 
of X using the formulas: 



x s = < 



i £ 1 < s < M, 

M 

it J2 %i,s-i+l, M<s<k, (8) 

1=1 

N-s+1 

TV-s+i XI ^i+s-k,k-i+lj k < S < N. 
i=l 

4. MATHEMATICAL MODELLING 

We computed a mathematical model to analyse the accuracy of the component recon- 
struction using the SSA method. We took the pulse-TOA residuals of the pulsars PSR 
J0613-0200, J1640+2224, J16431224, J1713+0747, J1939+2134, and J2145-0750 as the ini- 
tial reconstructed series. We added to these a series of normally distributed random values 
with zero mean and a dispersion based on the instrumental error in the pulse TOAs. In total, 
we generated 200 random realizations for each pulsar. The instrumental errors were taken 
from the paper of Oreshko [11], who presents the errors of the pulse TOAs determined using 
the AS-600 installation. The main input to the error is made by instability of the frequency 
of the bandpass filter. For the AS-600, this is given by the formula a t (f) = 0.0073 DM /is. 
Additional sources of inaccuracy are the local time standard (~ 100 ns) and the period 
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synthesizer ((~ 10 ns) and recorder (~ 20 ns). Table 2 gives the instrumental errors a T for 
the pulse TOAs for pulsars with various DMs and the mean error of each decomposition 
component <j\ . 

Next, the series of pulse TOAs were decomposed into components using the SSA method 
and the rms difference between the reconstructed and initial components was calculated. 
Figure 2 shows the results of our modelling for each of the six pulsars. The rms deviations 
are not uniformly distributed between the components - on average, components with higher 
numbers have lower errors. The error in the reconstructed components comprises 19%-32% 
of the instrumental error. 

We computed another mathematical model to analyse the real accuracy of the angular 
correlation function (1). Six pulsars were randomly placed on the celestial sphere in a 
Cartesian coordinate system Oxyz [2]. Next, the relative variation of the pulse frequency 
was calculated according to the relative positions of the pulsar and the gravitational wave, 
and the angular correlation coefficient ((9) was calculated. In all, N = 51 trials were made, in 
accordance with the number of data for each pulsar. Figure 3 shows the two-point correlation 
function obtained from these calculations averaged over 0.13 rad. The random error in the 
correlation coefficient is a p = 0.15. 

5. RESULTS 

To aid a visual analysis, we artificially added a seventh "pulsar" with the mean residuals 
[16] and the mean coordinates of the six pulsars. This does not add new data, but improves 
the filling of the plot. 

Figure 4 shows the component decomposition of the pulse-TOA residuals using the SSA 
method. The decrease in the period and amplitude of the components as the component 
number increases is clearly visible. We calculated power spectra (periodograms) for each 
component, shown in Fig. 5. An interesting feature of these spectra is that components 
with low numbers (1-5) display a well defined rise at low frequencies, corresponding to 
the presence of so-called "red" noises. Components 6-9 display a maximum in the middle 
of the frequency range, while components 10-12 rise at high frequencies (so-called "blue" 
spectra) . Thus, in terms of spectral indices, the spectral index of the periodogram grows as 
the component number increases. Since, generally speaking, spectra with variable spectral 
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indices can be associated with noise having various physical natures [5], we conclude that 
the SSA method for decomposing the time series into components enables us to distinguish 
components having various physical origins. 

We calculated the angular correlation functions given by (1) for the decomposition compo- 
nents of a given number, shown in Fig. 6. The experimental two-point correlation functions 
for components 8 and 9 shows a good agreement with the theoretical function: the correla- 
tion coefficient between the theoretical and experimental points is p = 0.92 ± 0.10. For this 
correlation coefficient and the number of points N = 19, the probability that this correlation 
occurred by chance is Pr(p = 0.92, N = 19) < 10~ 7 [17]. The two-point correlation function 
for the 8th component smoothed over two to five points is shown in more detail in Fig. 7. 

There are several possible explanations for this result. The most interesting would be 
the actual detection of gravitational waves. If the correlated amplitude of the 8th and 9th 
components of the decomposition, equal to 0.5 fis, is recalculated to the energy density 
of the inferred gravitational-wave background taking into account the observation span of 
six years [18], we obtain Q g h 2 ~ 10~ 8 , which is much higher than previous upper limits 
for the gravitational- wave background obtained in many studies (see, e.g., [19]), which give 
Q g h 2 < 10- 10 . 

It is also possible that this result is associated with some as yet unknown process, which 
behaves similarly to a gravitational- wave background with respect to the angular correlation 
function, and which requires additional study. 

6. CONCLUSION 

We have suggested a new method for detecting a gravitational-wave background based 
on observations of a group of millisecond pulsars. The method employs a new approach, 
based on a preliminary decomposition of the pulse TOA residuals of pulsars into spectral 
components with various physical natures and their subsequent analysis using the angular 
correlation function. We have applied this approach to the pulse-TOA residuals for the six 
pulsars PSR J0613-0200, J1640+2224, J1643-1224, J1713+0747, J1939+2134, and J2145- 
0750. This has yielded a calculated angular cross-correlation function that agrees with the 
function derived for the gravitational- wave background at a statistically significant level. 
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Table 1. Estimates of the timing parameters for the six pulsars 
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Table 2. Instrumental errors in the pulse TO As a T as functions of the DM and the error in the 
reconstruction of the decomposition components a\ . 
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Figure 1. Pulse-TOA residuals for the six pulsars averaged over 40 days (in /zs). The modified 
Julian date is plotted along the horizontal axis. 




Figure 2. Error in the reconstruction u\ (y axis) in fis as a function of the number of 
decomposition components k (x axis) for the six pulsars. 




Figure 3. Results of numerical modelling for the two-point correlation function. The data for the 

six pulsars with TV = 51 pulse TOAs is used. 
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Figure 4. Decomposition of TOA residuals into components for the six pulsars using the SSA 
method. The pulsar number and component number are indicated above each plot. The horizontal 
axis plots the time in 40-day intervals, and the vertical axis the time in fis. 
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Figure 5. Power spectra of the TOA-residual decomposition components for the six pulsars 
obtained using the SSA method. The pulsar number and component number are indicated above 
each plot. The horizontal axis plots the frequency (in yr _1 ) and the vertical axis the power (in /xs 2 ). 



Figure 6. Angular-correlation function calculated for the various TOA-residual decomposition 
components for the six pulsars smoothed over three points (the component numbers are indicated 
above the plots). The solid curves show the theoretical dependencies calculated assuming the 
presence of a gravitational-wave background. Components 8 and 9 show statistically significant 

correlations with the theoretical curves. 




Figure 7. Angular-correlation function calculated for the eighth component of the TOA-residual 
decomposition for the six pulsars smoothed over two to five points. The solid curve shows the 
theoretical dependence calculated assuming the presence of a gravitational-wave background. 



